Enteropathogen antibody dynamics and force of infection among children in low-resource settings
Kenya antibody profiles, Figure 4 - supplements 1 and 3
Notebook Summary
This script plots individual antibody trajectories in the Asembo, Kenya cohort from enrollment to follow-up. 205 children were measured at both timepoints.
Script preamble
## [1] "/Users/benarnold/enterics-seroepi"Identify incident changes
Identify children whose status changed between enrollment and follow-up. Those who changed from negative to positive are seroconverters (seroi below), and those who changed from positive to negative are seroreverters (seror below).
#-----------------------------
# identify incident
# seroconversions and reversions
# based on crossing the
# seropositivity cutoff
# and
# based on a 4-fold change in MFI
# to above the cutoff (seroconversion)
# or starting above cutoff (seroreversion)
#-----------------------------
dlp <- dl %>%
group_by(antigen,antigenf,childid) %>%
mutate(seroposA=ifelse(time=="A",seropos,NA),
seroposA=max(seroposA,na.rm=TRUE),
seroposB=ifelse(time=="B",seropos,NA),
seroposB=max(seroposB,na.rm=TRUE),
seroi=ifelse(seroposB==1 & seroposA==0,1,0),
seror=ifelse(seroposB==0 & seroposA==1,1,0),
mfiA=ifelse(time=="A",logmfi,NA),
mfiA=max(mfiA,na.rm=TRUE),
mfiB=ifelse(time=="B",logmfi,NA),
mfiB=max(mfiB,na.rm=TRUE),
seroi4fold=ifelse(logmfidiff>log10(4) & mfiB>serocut,1,0),
seror4fold=ifelse(logmfidiff< - log10(4) & mfiA>serocut,1,0)
) Antibody MFI trajectories
Longitudinal antibody trajectories of individual children for all antigens measured, colored by the change in antibody levels (increases and decreases)
#-----------------------------
# Plot individual antibody
# trajectories, colored
# by increases and decreases
#-----------------------------
# create an change category factor for plot aesthetics
dlp <- dlp %>%
mutate(mficat=factor(ifelse(logmfidiff>0,"Increase","Decrease"),levels=c("Increase","Decrease")))
# convert time to numeric (for plotting)
dlp <- dlp %>%
mutate(svy=ifelse(time=="A",0,1))
log10labs <- c(
expression(10^0),
expression(10^1),
expression(10^2),
expression(10^3),
expression(10^4)
)
ppq <- ggplot(data=filter(dlp,!is.na(mficat)),aes(x=svy,y=logmfi,group=factor(childid),color=mficat)) +
# geom_point(alpha=0.2) +
geom_line(alpha=0.2) +
# geom_hline(aes(yintercept=mixcut),linetype="dashed") +
scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
scale_y_continuous(breaks=0:4,labels=log10labs) +
scale_color_manual(values=c(corange,cteal,"gray80"),
guide=guide_legend(title="MFI change:",override.aes = list(alpha=1))
) +
facet_wrap(~antigenf,nrow=6,ncol=2) +
labs(x="Measurement",y="Luminex response (MFI-bg)") +
theme_minimal() +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
axis.ticks.y=element_line(color="gray40"),
legend.position = "top"
)
ppq
Create the same figure of individual antibody trajectories between measurements, but limit it to antigens with seroincidence measures, and color the trajectories by seroconversion/reversion status.
#-----------------------------
# pull the incidence info
# and just merge it back
# to the full longidudinal data
#-----------------------------
dl2 <- dlp
# create an incidence category factor for plot aesthetics
# for seroconversion based on seropositivity only
dl2$icat <- factor(NA,levels=c("Seroconversion","Seroreversion","No change"))
dl2$icat[dl2$seroi==1] <- "Seroconversion"
dl2$icat[dl2$seror==1] <- "Seroreversion"
dl2$icat[dl2$seroi==0 & dl2$seror==0] <- "No change"
# create an incidence category factor for plot aesthetics
# for seroconversion based on 4-fold changes in MFI
dl2$i4cat <- factor(NA,levels=c("Seroconversion","Seroreversion","No change"))
dl2$i4cat[dl2$seroi4fold==1] <- "Seroconversion"
dl2$i4cat[dl2$seror4fold==1] <- "Seroreversion"
dl2$i4cat[dl2$seroi4fold==0 & dl2$seror==0] <- "No change"
# create an incidence category factor for plot aesthetics
# for seroconversion comparing both methods
dl2$icat_comp <- factor(NA,levels=c(">4-fold increase, across cutoff",">4-fold decrease",">4-fold increase, above cutoff","<4-fold change"))
dl2$icat_comp[dl2$seroi4fold==1 & dl2$seroi==1] <- ">4-fold increase, across cutoff"
dl2$icat_comp[dl2$seror4fold==1] <- ">4-fold decrease"
dl2$icat_comp[dl2$seroi4fold==1 & dl2$seroi==0] <- ">4-fold increase, above cutoff"
dl2$icat_comp[is.na(dl2$icat_comp)] <- "<4-fold change"
dl2 <- dl2 %>%
mutate(icat_comp = factor(icat_comp,levels=c(">4-fold increase, across cutoff",
">4-fold increase, above cutoff",
">4-fold decrease",
"<4-fold change"))
)
pp <- ggplot(data=filter(dl2,!is.na(icat)),aes(x=svy,y=logmfi,group=factor(childid),color=icat_comp)) +
# geom_point(alpha=0.2) +
geom_hline(aes(yintercept=serocut),linetype="dashed",size=0.3) +
geom_line(alpha=0.2) +
scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
scale_y_continuous(breaks=0:4,labels=log10labs) +
scale_color_manual(values=c(corange,cmagent,cteal,"gray80"),
guide=guide_legend(title="MFI change:", override.aes = list(alpha=1), ncol=2,nrow=2),
) +
facet_wrap(~antigenf,nrow=6,ncol=2) +
labs(x="Measurement",y="Luminex response (MFI-bg)") +
theme_minimal() +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
axis.ticks.y=element_line(color="gray40"),
legend.position = "top"
)
pp
Figure 4 - supplement 1
#-----------------------------
# final figure, excluding cholera
# due to cross-reactivity with ETEC
#-----------------------------
relev <- levels(dl2$antigenf)[c(1:6,11,7,9,10)]
dl3 <- dl2 %>%
filter(!is.na(icat) & antigen!="cholera") %>%
ungroup() %>%
mutate(antigenf=factor(antigenf,levels=relev))
pp2 <- ggplot(data=dl3,aes(x=svy,y=logmfi,group=factor(childid),color=icat_comp)) +
# geom_point(alpha=0.2) +
geom_hline(aes(yintercept=serocut),linetype="dashed",size=0.3) +
geom_line(alpha=0.2) +
scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
scale_y_continuous(breaks=0:4,labels=log10labs) +
scale_color_manual(values=c(corange,cteal,cmagent,"gray80"),
guide=guide_legend(title="MFI change:", override.aes = list(alpha=1), ncol=2,nrow=2),
) +
facet_wrap(~antigenf,nrow=6,ncol=2) +
labs(x="Measurement",y="Luminex response (MFI-bg)") +
theme_minimal() +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
axis.ticks.y=element_line(color="gray40"),
legend.position = "top"
)
pp2
Figure 4 - supplement 3
Estimate prevalence of each category
dl4 <- dl3 %>%
group_by(antigenf,icat_comp, .drop=FALSE) %>%
tally() %>%
group_by(antigenf) %>%
mutate(N=sum(n),
q=N-n,
prev=n/N
)
# exact binomial confidence intervals on prevalence estimates
ci_lb <- apply(dl4[c("n","q")],1,function(x) binom.test(x,alternative="two.sided")$conf.int[1])
ci_ub <- apply(dl4[c("n","q")],1,function(x) binom.test(x,alternative="two.sided")$conf.int[2])
dl4 <- dl4 %>%
bind_cols(data.frame(ci_lb,ci_ub)) %>%
arrange(antigenf,icat_comp)# compare seroconversion classification prevalences by measurement
sero_comp_p <- ggplot(data=dl4, aes(x=icat_comp, y=prev,color=icat_comp) ) +
facet_wrap(~antigenf,nrow=6,ncol=2)+
geom_errorbar(aes(ymin=ci_lb,ymax=ci_ub),width=0.1)+
geom_point()+
scale_color_manual(values=c(corange,cmagent,cteal,"gray50"),
guide=guide_legend(title="MFI change:",
override.aes = list(alpha=1),
nrow=2,ncol=2
))+
xlab("From Enrollment to Follow-up") +
ylab("Percentage of children (%)") +
coord_cartesian(ylim=c(0,1),xlim=c(1,4))+
scale_y_continuous(breaks=seq(0,1,by=0.2),labels=sprintf("%1.0f",seq(0,1,by=0.2)*100)) +
# scale_x_continuous(breaks=1:5) +
theme_minimal() +
theme(
legend.position="top",
panel.grid.minor.x=element_blank(),
axis.text.x = element_blank()
)
sero_comp_p
Session Info
## R version 3.5.3 (2019-03-11)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ROCR_1.0-7 gplots_3.0.1 tmleAb_0.3.4
## [4] nlme_3.1-137 gridExtra_2.3 RColorBrewer_1.1-2
## [7] ellipse_0.4.1 scales_1.0.0 knitr_1.22
## [10] mixtools_1.1.0 viridis_0.5.1 viridisLite_0.3.0
## [13] foreign_0.8-71 lubridate_1.7.4 doParallel_1.0.11
## [16] iterators_1.0.9 foreach_1.4.4 kableExtra_0.8.0.0001
## [19] forcats_0.3.0 stringr_1.4.0 dplyr_0.8.0.1
## [22] purrr_0.3.2 readr_1.1.1 tidyr_0.8.3
## [25] tibble_2.1.1 ggplot2_3.1.1 tidyverse_1.2.1
## [28] here_0.1
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.0 jsonlite_1.6 splines_3.5.3
## [4] modelr_0.1.2 gtools_3.5.0 assertthat_0.2.1
## [7] highr_0.8 cellranger_1.1.0 yaml_2.2.0
## [10] pillar_1.4.0 backports_1.1.4 lattice_0.20-38
## [13] glue_1.3.1 digest_0.6.18 rvest_0.3.2
## [16] colorspace_1.3-2 htmltools_0.3.6 Matrix_1.2-15
## [19] plyr_1.8.4 psych_1.8.4 pkgconfig_2.0.2
## [22] broom_0.4.4 haven_2.1.0 gdata_2.18.0
## [25] withr_2.1.2 lazyeval_0.2.2 cli_1.1.0
## [28] mnormt_1.5-5 survival_2.43-3 magrittr_1.5
## [31] crayon_1.3.4 readxl_1.1.0 evaluate_0.13
## [34] MASS_7.3-51.1 segmented_0.5-3.0 xml2_1.2.0
## [37] tools_3.5.3 hms_0.4.2 munsell_0.5.0
## [40] packrat_0.4.9-3 compiler_3.5.3 caTools_1.17.1
## [43] rlang_0.3.4 grid_3.5.3 rstudioapi_0.9.0
## [46] bitops_1.0-6 labeling_0.3 rmarkdown_1.12
## [49] gtable_0.3.0 codetools_0.2-16 reshape2_1.4.3
## [52] R6_2.4.0 rprojroot_1.3-2 KernSmooth_2.23-15
## [55] stringi_1.4.3 Rcpp_1.0.1 tidyselect_0.2.5
## [58] xfun_0.6